// windows related: to show progress status
#include "mcmcGibbsLearning.h"


//extern CProgressCtrl *pCtrProgressActive;

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
// Maximum psedo-likelihood estimator

// constructoions and destructions
MPLE_Iid::MPLE_Iid()
{
}

MPLE_Iid::~MPLE_Iid()
{
}

// operations
void MPLE_Iid::Set(const int num_label1,const int nlabel1, mcmcMatrix<int> &mx_label,const int num_f, const irmm_Rect bd)
{
	num_label	= num_label1;
	nlabel		= nlabel1;
	pmx_label	= &mx_label;
	bound		= bd;
	num_filters = num_f;

	mx_beta.SetDimension(num_label,num_filters);
	mx_histogram.SetDimension(num_label,num_filters);
	mx_histogram_syn.SetDimension(num_label,num_filters);
	mx_Z.SetDimension(num_label,num_filters);
}

void MPLE_Iid::SetEachFilter(const int ifilter, Raster<double> &ras, double min_b, double max_b, double delta)
{
	pras[ifilter] = &ras;
	for (int i=0; i<num_label; i++)
	{
		mx_beta(i,ifilter).Set(min_b,max_b,delta);
		mx_histogram(i,ifilter).Set(min_b,max_b,delta);
		mx_histogram(i,ifilter).Init(0.0);
		mx_histogram_syn(i,ifilter).Set(min_b,max_b,delta);
	}
}

void MPLE_Iid::Learn(const int time)
{
	int		i,j,k,m,t;
	mcmcVector<double> vc_temp;

	// compute histogram first
	// compute the mean value
	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && (*pmx_label)(j,i)>=nlabel && (*pmx_label)(j,i)<=nlabel+num_label-1)
			{
				k = (*pmx_label)(j,i)-nlabel;

				for (m=0; m<num_filters; m++)
					mx_histogram(k,m).AddBin(pras[m]->Data(j,i),1.0);
			}
		}
	// normalization
	for (k=0; k<num_label; k++)
		for (m=0; m<num_filters; m++)
			mx_histogram(k,m).Normalize();

	// initialization
	for (k=0; k<num_label; k++)
		for (m=0; m<num_filters; m++)
		{
			mx_beta(k,m).Init(1.0);
			mx_Z(k,m) = mx_beta(k,m).ComputeZ();
		}

	double d;

	for (t=0; t<time; t++)
	{
		if (t<300)			d = 0.1;
		else if (t<800)		d = 0.05;
		else				d = 0.01;

		for (k=0; k<num_label; k++)
			for (m=0; m<num_filters; m++)
			{
				mx_histogram_syn(k,m).Init(0.0);
				for (i=0; i<mx_histogram(k,m).count; i++)
					mx_histogram_syn(k,m).pbins[i] = 1.0/mx_Z(k,m)*exp(-mx_beta(k,m).pbins[i]);
				mx_histogram_syn(k,m).Normalize();

				for (i=0; i<mx_beta(k,m).count; i++)
					mx_beta(k,m).pbins[i] = mx_beta(k,m).pbins[i]-(mx_histogram(k,m).pbins[i]-mx_histogram_syn(k,m).pbins[i])*d;
				// compute Z again
				mx_Z(k,m) = mx_beta(k,m).ComputeZ();
			}
	}
}


// construction & destruction
MPLE_Iid_Luv::MPLE_Iid_Luv()
{
}

MPLE_Iid_Luv::~MPLE_Iid_Luv()
{
}

void MPLE_Iid_Luv::Set(const int num_label1,const int nlabel1, mcmcMatrix<int> &mx_label,const irmm_Rect bd)
{
	num_label	= num_label1;
	nlabel		= nlabel1;
	pmx_label	= &mx_label;
	bound		= bd;
	num_filters = 6;

	mx_beta.SetDimension(num_label,num_filters);
	mx_histogram.SetDimension(num_label,num_filters);
	mx_histogram_syn.SetDimension(num_label,num_filters);
	mx_Z.SetDimension(num_label,1);
}

void MPLE_Iid_Luv::SetData(Raster<double> &rasL, Raster<double> &rasU, Raster<double> &rasV,
						   Raster<double> &rasLU, Raster<double> &rasLV, Raster<double> &rasUV,
						   mcmcVector<double> &vc_min_b1, mcmcVector<double> &vc_max_b1, mcmcVector<double> &vc_delta1)
{
	pras[0]		= &rasL;
	pras[1]		= &rasU;
	pras[2]		= &rasV;
	pras[3]		= &rasLU;
	pras[4]		= &rasLV;
	pras[5]		= &rasUV;

	vc_min_b	= vc_min_b1;
	vc_max_b	= vc_max_b1;
	vc_delta	= vc_delta1;

	for (int i=0; i<num_label; i++)
	{
		for (int j=0; j<num_filters; j++)
		{
			mx_beta(i,j).Set(vc_min_b(j),vc_max_b(j),vc_delta(j));
			mx_histogram(i,j).Set(vc_min_b(j),vc_max_b(j),vc_delta(j));
			mx_histogram(i,j).Init(0.0);
			mx_histogram_syn(i,j).Set(vc_min_b(j),vc_max_b(j),vc_delta(j));
		}
	}
}

void MPLE_Iid_Luv::Learn(const int time)
{
	int				i,j,k,m,t;
	double			l,u,v;
	double			d,prob;
	mcmcVector<double>	vc_temp;

	// compute histogram first
	// compute the mean value
	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if ((*pmx_label)(j,i)>0 && (*pmx_label)(j,i)>=nlabel && (*pmx_label)(j,i)<=nlabel+num_label-1)
			{
				k = (*pmx_label)(j,i)-nlabel;

				for (m=0; m<num_filters; m++)
					mx_histogram(k,m).AddBin(pras[m]->Data(j,i),1.0);
			}
		}
	// normalization
	for (k=0; k<num_label; k++)
		for (m=0; m<num_filters; m++)
		{
			mx_histogram(k,m).Normalize();
			mx_beta(k,m).Init(1.0);
		}

	for (t=0; t<time; t++)
	{
		if (t<300)			d = 0.1;
		else if (t<800)		d = 0.05;
		else				d = 0.01;

		for (k=0; k<num_label; k++)
		{
			for (m=0; m<num_filters; m++)
				mx_histogram_syn(k,m).Init(0.0);

			for (l=vc_min_b(0); l<=vc_max_b(0); l=l+vc_delta(0))
			for (u=vc_min_b(1); u<=vc_max_b(1); u=u+vc_delta(1))
			for (v=vc_min_b(2); v<=vc_max_b(2); v=v+vc_delta(2))
			{
				prob = Prob(k,l,u,v,false);

				mx_histogram_syn(k,0).AddBin(l,prob);
				mx_histogram_syn(k,1).AddBin(u,prob);
				mx_histogram_syn(k,2).AddBin(v,prob);
				mx_histogram_syn(k,3).AddBin(l+u,prob);
				mx_histogram_syn(k,4).AddBin(l+v,prob);
				mx_histogram_syn(k,5).AddBin(u+v,prob);
			}

			for (m=0; m<num_filters; m++)
			{
				mx_histogram_syn(k,m).Normalize();
				for (int ii=0; ii<mx_beta(k,m).count; ii++)
					mx_beta(k,m).pbins[ii] = mx_beta(k,m).pbins[ii]
						-(mx_histogram(k,m).pbins[ii]-mx_histogram_syn(k,m).pbins[ii])*d;
			}

		}
	}

	mx_Z.SetDimension(num_label,1);
	mx_Z.InitValue(0.0);
	// compute Z
	for (k=0; k<num_label; k++)
	{
		for (l=vc_min_b(0); l<=vc_max_b(0); l=l+vc_delta(0))
		for (u=vc_min_b(1); u<=vc_max_b(1); u=u+vc_delta(1))
		for (v=vc_min_b(2); v<=vc_max_b(2); v=v+vc_delta(2))
			mx_Z(k) += Prob(k,l,u,v,false);

	}
}


void MPLE_Iid_Luv::Synthesiz(Raster_Color<uchar> &image_color)
{
	RF3D	rf_sampler;
	double	pLUV[3];
	int		j,i;
	double	l,u,v;

	rf_sampler.Set(vc_min_b(0),vc_max_b(0),vc_delta(0),
				   vc_min_b(1),vc_max_b(1),vc_delta(1),
				   vc_min_b(2),vc_max_b(2),vc_delta(2));

	for (l=vc_min_b(0); l<=vc_max_b(0); l=l+vc_delta(0))
	for (u=vc_min_b(1); u<=vc_max_b(1); u=u+vc_delta(1))
	for (v=vc_min_b(2); v<=vc_max_b(2); v=v+vc_delta(2))
		rf_sampler.SetOneBin(l,u,v,Prob(0,l,u,v,true));
	rf_sampler.SetCdf();

	for (j=0; j<image_color.rows(); j++)
		for (i=0; i<image_color.cols(); i++)
		{
			rf_sampler.Sampling(pLUV[0],pLUV[1],pLUV[2]);
			image_color.SetPixelInLUV(j,i,(int)pLUV[0],(int)pLUV[1],(int)pLUV[2]);
		}
}

/////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////
MPLE_Luv::MPLE_Luv()
{
	bprecompute_resp = true;
	for (int i=0; i<6; i++)
		p_important_color_index[i] = NULL;
	for (int i=0; i<3; i++)
		p_important_color[i] = NULL;
}

MPLE_Luv::~MPLE_Luv()
{
	for (int i=0; i<6; i++)
		if (p_important_color_index[i] != NULL)
			delete [](p_important_color_index[i]);
	for (int i=0; i<3; i++)
		if (p_important_color[i] != NULL)
			delete [](p_important_color[i]);
}

bool MPLE_Luv::ValidLUV(const double l, const double u, const double v)
{
	int k_l,k_u,k_v,k_lu,k_lv,k_uv;

	k_l  = frame_channel_l.GetIndex(l);
	k_u  = frame_channel_u.GetIndex(u);
	k_v  = frame_channel_v.GetIndex(v);
	k_lu = frame_channel_lu.GetIndex(l+u);
	k_lv = frame_channel_lv.GetIndex(l+v);
	k_uv = frame_channel_uv.GetIndex(u+v);

	if (frame_channel_l.ValidIndex(k_l) && frame_channel_u.ValidIndex(k_u) &&
		frame_channel_v.ValidIndex(k_v) && frame_channel_lu.ValidIndex(k_lu) &&
		frame_channel_lv.ValidIndex(k_lv) && frame_channel_uv.ValidIndex(k_uv))
		return true;
	else
		return false;
}

void MPLE_Luv::Initialization(void)
{
	// learn color texture
	// L
	frame_channel_l.Set(presp->filter_bank_L,true,biggest_rowsize,biggest_colsize);
	// U
	frame_channel_u.Set(presp->filter_bank_U,true,biggest_rowsize,biggest_colsize);
	// V
	frame_channel_v.Set(presp->filter_bank_V,true,biggest_rowsize,biggest_colsize);

	// get all the existing colors as important ones for computing pseudo likelihood probability
	int		i,j,k,m;
	double	l,u,v;
	Histogram *phist_l,*phist_u,*phist_v;

	phist_l = &frame_channel_l.histogram_ref_input;
	phist_u = &frame_channel_u.histogram_ref_input;
	phist_v = &frame_channel_v.histogram_ref_input;

	num_total_color = 0;
	for (j=0; j<phist_l->count; j++)
	{
		l = phist_l->GetValueIndex(j);
		for (i=0; i<phist_u->count; i++)
		{
			u = phist_u->GetValueIndex(i);
			for (k=0; k<phist_v->count; k++)
			{
				v = phist_v->GetValueIndex(k);
				if (ValidLUV(l,u,v))
					num_total_color++;
			}
		}
	}

	// L,U,V are stored repsectively
	m = 0;
	for (i=0; i<6; i++)
		p_important_color_index[i] = new int[num_total_color];
	for (i=0; i<3; i++)
		p_important_color[i] = new double[num_total_color];

	for (j=0; j<phist_l->count; j++)
	{
		l = phist_l->GetValueIndex(j);
		for (i=0; i<phist_u->count; i++)
		{
			u = phist_u->GetValueIndex(i);
			for (k=0; k<phist_v->count; k++)
			{
				v = phist_v->GetValueIndex(k);
				if (ValidLUV(l,u,v))
				{
					p_important_color[0][m] = l;
					p_important_color[1][m] = u;
					p_important_color[2][m] = v;

					p_important_color_index[0][m] = frame_channel_l.GetIndex(l);
					p_important_color_index[1][m] = frame_channel_u.GetIndex(u);
					p_important_color_index[2][m] = frame_channel_v.GetIndex(v);
					p_important_color_index[3][m] = frame_channel_lu.GetIndex(l+u);
					p_important_color_index[4][m] = frame_channel_lv.GetIndex(l+v);
					p_important_color_index[5][m] = frame_channel_uv.GetIndex(u+v);

					m++;
				}
			}
		}
	}
}

bool MPLE_Luv::UpdateBeta(const double d)
{
	bool bConverge=true;

	// L
	bConverge = frame_channel_l.UpdateBeta(d) & bConverge;
	// U
	bConverge = frame_channel_u.UpdateBeta(d) & bConverge;
	// V
	bConverge = frame_channel_v.UpdateBeta(d) & bConverge;
	// L+U
	bConverge = frame_channel_lu.UpdateBeta(d) & bConverge;
	// L+V
	bConverge = frame_channel_lv.UpdateBeta(d) & bConverge;
	// U+V
	bConverge = frame_channel_uv.UpdateBeta(d) & bConverge;

	return bConverge;
}

void MPLE_Luv::ComputeSiteHistogramInAdvance(const int kcluster)
{
	if (bprecompute_resp)
	{
		// compute histogram in advance
		frame_channel_l.ComputeResponse(presp->cut_off,presp->rows()-(presp->cut_off+1),
										presp->cut_off,presp->cols()-(presp->cut_off+1));
		frame_channel_u.ComputeResponse(presp->cut_off,presp->rows()-(presp->cut_off+1),
										presp->cut_off,presp->cols()-(presp->cut_off+1));
		frame_channel_v.ComputeResponse(presp->cut_off,presp->rows()-(presp->cut_off+1),
										presp->cut_off,presp->cols()-(presp->cut_off+1));
		frame_channel_lu.ComputeResponse(presp->cut_off,presp->rows()-(presp->cut_off+1),
										presp->cut_off,presp->cols()-(presp->cut_off+1));
		frame_channel_lv.ComputeResponse(presp->cut_off,presp->rows()-(presp->cut_off+1),
										presp->cut_off,presp->cols()-(presp->cut_off+1));
		frame_channel_uv.ComputeResponse(presp->cut_off,presp->rows()-(presp->cut_off+1),
										presp->cut_off,presp->cols()-(presp->cut_off+1));

		SaveRespHistToFile(0);
	}
}

void MPLE_Luv::SaveRespHistToFile(const int m)
{
	char		str[200];
	ofstream	file;

	sprintf(str,"d:\\temp\\%d.rsp", m);

	file.open(str, ios::out|ios::binary);
	if (file.fail())
		return;

	frame_channel_l.SaveRespHistToFile(file);
	frame_channel_u.SaveRespHistToFile(file);
	frame_channel_v.SaveRespHistToFile(file);
	frame_channel_lu.SaveRespHistToFile(file);
	frame_channel_lv.SaveRespHistToFile(file);
	frame_channel_uv.SaveRespHistToFile(file);

	file.close();
}

void MPLE_Luv::ReadRespHistFromFile(const int m)
{
	if (current_patch != m)
	{
		char	str[200];
		ifstream	file;

		sprintf(str,"d:\\temp\\%d.rsp", m);

		file.open(str, ios::in|ios::binary|ios::_Nocreate);
		if (file.fail())
			return;

		frame_channel_l.ReadRespHistFromFile(file);
		frame_channel_u.ReadRespHistFromFile(file);
		frame_channel_v.ReadRespHistFromFile(file);
		frame_channel_lu.ReadRespHistFromFile(file);
		frame_channel_lv.ReadRespHistFromFile(file);
		frame_channel_uv.ReadRespHistFromFile(file);

		file.close();

		current_patch = m;
	}
}

double	MPLE_Luv::Energy(const double l, const double u, const double v)
{
	int			k;
	double		energy;

	energy = 0.0;
	// L
	k = frame_channel_l.GetIndex(l);
	if (frame_channel_l.vc_energy.Valid(k))
		energy += frame_channel_l.vc_energy(k);
	else
		return 9000.0;

	// U
	k = frame_channel_u.GetIndex(u);
	if (frame_channel_u.vc_energy.Valid(k))
		energy += frame_channel_u.vc_energy(k);
	else
		return 9000.0;
	// V
	k = frame_channel_v.GetIndex(v);
	if (frame_channel_v.vc_energy.Valid(k))
		energy += frame_channel_v.vc_energy(k);
	else
		return 9000.0;
	// L+U
	k = frame_channel_lu.GetIndex(l+u);
	if (frame_channel_lu.vc_energy.Valid(k))
		energy += frame_channel_lu.vc_energy(k);
	else
		return 9000.0;
	// L+V
	k = frame_channel_lv.GetIndex(l+v);
	if (frame_channel_lv.vc_energy.Valid(k))
		energy += frame_channel_lv.vc_energy(k);
	else
		return 9000.0;
	// U+V
	k = frame_channel_uv.GetIndex(u+v);
	if (frame_channel_uv.vc_energy.Valid(k))
		energy += frame_channel_uv.vc_energy(k);
	else
		return 9000.0;

	return energy;
}

double MPLE_Luv::Energy(const int k_l, const int k_u, const int k_v, const int k_lu,const int k_lv,const int k_uv)
{
	double		energy;

	energy = 0.0;
	// L
	energy += frame_channel_l.vc_energy(k_l);
	// U
	energy += frame_channel_u.vc_energy(k_u);
	// V
	energy += frame_channel_v.vc_energy(k_v);
	// L+U
	energy += frame_channel_lu.vc_energy(k_lu);
	// L+V
	energy += frame_channel_lv.vc_energy(k_lv);
	// U+V
	energy += frame_channel_uv.vc_energy(k_uv);

	return energy;
}


void MPLE_Luv::UpdateHistogramSyn(void)
{
	// L
	frame_channel_l.UpdateHistogramSyn();
	// U
	frame_channel_u.UpdateHistogramSyn();
	// V
	frame_channel_v.UpdateHistogramSyn();
	// L+U
	frame_channel_lu.UpdateHistogramSyn();
	// L+V
	frame_channel_lv.UpdateHistogramSyn();
	// U+V
	frame_channel_uv.UpdateHistogramSyn();
}

// from precomputed histogram
void MPLE_Luv::UpdateHistogramSynTemp(const double prob, const double l, const double u, const double v, const int row, const int col)
{
	// L
	frame_channel_l.UpdateHistogramSynTemp(prob,l,row,col);
	// U
	frame_channel_u.UpdateHistogramSynTemp(prob,u,row,col);
	// V
	frame_channel_v.UpdateHistogramSynTemp(prob,v,row,col);
	// L+U
	frame_channel_lu.UpdateHistogramSynTemp(prob,l+u,row,col);
	// L+V
	frame_channel_lv.UpdateHistogramSynTemp(prob,l+v,row,col);
	// U+V
	frame_channel_uv.UpdateHistogramSynTemp(prob,u+v,row,col);
}

void MPLE_Luv::ClearHistogramSyn(void)
{
	// L
	frame_channel_l.ClearHistogramSyn();
	// U
	frame_channel_u.ClearHistogramSyn();
	// V
	frame_channel_v.ClearHistogramSyn();
	// L+U
	frame_channel_lu.ClearHistogramSyn();
	// L+V
	frame_channel_lv.ClearHistogramSyn();
	// U+V
	frame_channel_uv.ClearHistogramSyn();
}

void MPLE_Luv::ClearHistogramSynTemp(void)
{
	// L
	frame_channel_l.ClearHistogramSynTemp();
	// U
	frame_channel_u.ClearHistogramSynTemp();
	// V
	frame_channel_v.ClearHistogramSynTemp();
	// L+U
	frame_channel_lu.ClearHistogramSynTemp();
	// L+V
	frame_channel_lv.ClearHistogramSynTemp();
	// U+V
	frame_channel_uv.ClearHistogramSynTemp();
}

void MPLE_Luv::UpdateFilterRespAtOneSite(const int j, const int i, const double l, const double u, const double v)
{
	// L
	frame_channel_l.UpdateFilterRespAtOneSite(j,i,l);
	// U
	frame_channel_u.UpdateFilterRespAtOneSite(j,i,u);
	// V
	frame_channel_v.UpdateFilterRespAtOneSite(j,i,v);
	// L+U
	frame_channel_lu.UpdateFilterRespAtOneSite(j,i,l+u);
	// L+V
	frame_channel_lv.UpdateFilterRespAtOneSite(j,i,l+v);
	// U+V
	frame_channel_uv.UpdateFilterRespAtOneSite(j,i,u+v);
}

void MPLE_Luv::Learn(const int time, const bool bprecompute,const int kcluster)
{
	int					i,j,m,t,win_size=20;
	register double		l,u,v;
	double				prob,energy,energy_min;
	mcmcVector<double>		vc_energy;

	// precompute all the responses regarding to the changes of pixel values
	bprecompute_resp = bprecompute;

	Initialization();
	// temporary energy
	vc_energy.SetRowVector(num_total_color);

	current_patch = -1;

	char btafile_name[200];
	sprintf(btafile_name,"a_%d.bta",kcluster);

	if (time > 0)
	{
		ComputeSiteHistogramInAdvance(kcluster);
		/* Jiadi
		if (pCtrProgressActive != NULL)
		{
			pCtrProgressActive->SetRange(0,time);
			pCtrProgressActive->SetPos(0);
			pCtrProgressActive->ShowWindow(SW_SHOW);
		}
		*/
	}

	register int mm;
	bool bConverge=false;
	t = 0;
	while (t<time && !bConverge)
	{
			for (m=0; m<10; m++)
			{
				ClearHistogramSyn();

				for (int ii=0; ii<win_size*win_size; ii++)
				{
					j = 0;
					i = 0;

					ClearHistogramSynTemp();
					energy_min = 1e10;

					// L
					frame_channel_l.SaveTempEnergy(j,i);	// from the pre-computed histog
					// U
					frame_channel_u.SaveTempEnergy(j,i);
					// V
					frame_channel_v.SaveTempEnergy(j,i);
					// L+U
					frame_channel_lu.SaveTempEnergy(j,i);
					// L+V
					frame_channel_lv.SaveTempEnergy(j,i);
					// U+V
					frame_channel_uv.SaveTempEnergy(j,i);

					// get the minimal energy value
					for (mm=0; mm<num_total_color; mm++)
					{
						energy = Energy(p_important_color_index[0][mm],p_important_color_index[1][mm],p_important_color_index[2][mm],
										p_important_color_index[3][mm],p_important_color_index[4][mm],p_important_color_index[5][mm]);
						vc_energy(mm) = energy;
							if (energy<energy_min)	energy_min = energy;
					}
					for (mm=0; mm<num_total_color; mm++)
					{
						l = p_important_color[0][mm];
						u = p_important_color[1][mm];
						v = p_important_color[2][mm];
							energy = vc_energy(mm)-energy_min;
						prob = ProbFromEnergy(energy);
							UpdateHistogramSynTemp(prob,l,u,v,j,i);
					}
				}

				UpdateHistogramSyn();
				// update betas
				bConverge = UpdateBeta(GetD(t));
				t++;
				if (t>=time)
					break;
				/* Jiadi
				if (pCtrProgressActive != NULL)
				{
					// show progress
					pCtrProgressActive->SetPos(t);
					pCtrProgressActive->ShowWindow(SW_SHOW);
				}
				*/
				if (bConverge)
					break;
			SaveBetaToFile(btafile_name);
			if (bConverge)
				break;
			}
	}
	/* Jiadi
	if (pCtrProgressActive != NULL)
		pCtrProgressActive->ShowWindow(SW_HIDE);
	*/

	if (time != 0)
		SaveBetaToFile(btafile_name);
	ReadBetaFromFile(btafile_name);
}

void MPLE_Luv::SaveBetaToFile(char *file_name)
{
	FILE		*pfile;
	Histogram	*phist;
	int			i;

	if (file_name == NULL)
		return;

	pfile = fopen(file_name,"w+");
	if(pfile != NULL)
	{
		for (i=0; i<frame_channel_l.vc_B.size(); i++)
		{
			phist = (Histogram *)&frame_channel_l.vc_B(i);
			phist->WriteToTxtFile(*pfile);
		}
		fprintf(pfile,"-----------------------\n");
		for (i=0; i<frame_channel_u.vc_B.size(); i++)
		{
			phist = (Histogram *)&frame_channel_u.vc_B(i);
			phist->WriteToTxtFile(*pfile);
		}
		fprintf(pfile,"-----------------------\n");
		for (i=0; i<frame_channel_v.vc_B.size(); i++)
		{
			phist = (Histogram *)&frame_channel_v.vc_B(i);
			phist->WriteToTxtFile(*pfile);
		}
		fprintf(pfile,"-----------------------\n");
		for (i=0; i<frame_channel_lu.vc_B.size(); i++)
		{
			phist = (Histogram *)&frame_channel_lu.vc_B(i);
			phist->WriteToTxtFile(*pfile);
		}
		fprintf(pfile,"-----------------------\n");
		for (i=0; i<frame_channel_lv.vc_B.size(); i++)
		{
			phist = (Histogram *)&frame_channel_lv.vc_B(i);
			phist->WriteToTxtFile(*pfile);
		}
		fprintf(pfile,"-----------------------\n");
		for (i=0; i<frame_channel_uv.vc_B.size(); i++)
		{
			phist = (Histogram *)&frame_channel_uv.vc_B(i);
			phist->WriteToTxtFile(*pfile);
		}
		fclose(pfile);
	}
}

void MPLE_Luv::ReadBetaFromFile(char *file_name)
{
	FILE		*pfile;
	Histogram	*phist;
	int			i;
	char		str[100];
	if (file_name == NULL)
		return;

	pfile = fopen(file_name,"r");
	if (pfile != NULL)
	{
		for (i=0; i<frame_channel_l.vc_B.size(); i++)
		{
			phist = (Histogram *)&frame_channel_l.vc_B(i);
			phist->ReadFromTxtFile(*pfile);
		}
		fscanf(pfile,"%s",&str);
		for (i=0; i<frame_channel_u.vc_B.size(); i++)
		{
			phist = (Histogram *)&frame_channel_u.vc_B(i);
			phist->ReadFromTxtFile(*pfile);
		}
		fscanf(pfile,"%s",&str);
		for (i=0; i<frame_channel_v.vc_B.size(); i++)
		{
			phist = (Histogram *)&frame_channel_v.vc_B(i);
			phist->ReadFromTxtFile(*pfile);
		}
		fscanf(pfile,"%s",&str);
		for (i=0; i<frame_channel_lu.vc_B.size(); i++)
		{
			phist = (Histogram *)&frame_channel_lu.vc_B(i);
			phist->ReadFromTxtFile(*pfile);
		}
		fscanf(pfile,"%s",&str);
		for (i=0; i<frame_channel_lv.vc_B.size(); i++)
		{
			phist = (Histogram *)&frame_channel_lv.vc_B(i);
			phist->ReadFromTxtFile(*pfile);
		}
		fscanf(pfile,"%s",&str);
		for (i=0; i<frame_channel_uv.vc_B.size(); i++)
		{
			phist = (Histogram *)&frame_channel_uv.vc_B(i);
			phist->ReadFromTxtFile(*pfile);
		}
		fclose(pfile);
	}
}

void MPLE_Luv::FilteringOnNewImage(Raster_Color<uchar> &image_color)
{
	Raster<double> ras_L,ras_U,ras_V,ras_LU,ras_LV,ras_UV;

	ras_L.SetDimension(image_color.rows(),image_color.cols());
	ras_U.SetDimension(image_color.rows(),image_color.cols());
	ras_V.SetDimension(image_color.rows(),image_color.cols());
	ras_LU.SetDimension(image_color.rows(),image_color.cols());
	ras_LV.SetDimension(image_color.rows(),image_color.cols());
	ras_UV.SetDimension(image_color.rows(),image_color.cols());

	for (int j=0; j<ras_L.rows(); j++)
		for (int i=0; i<ras_L.cols(); i++)
		{
			ras_L(j,i) = image_color.L(j,i);
			ras_U(j,i) = image_color.U(j,i);
			ras_V(j,i) = image_color.V(j,i);
			ras_LU(j,i) = image_color.L(j,i)+image_color.U(j,i);
			ras_LV(j,i) = image_color.L(j,i)+image_color.V(j,i);
			ras_UV(j,i) = image_color.U(j,i)+image_color.V(j,i);
		}

	// L
	frame_channel_l.FilteringOnNewImage(ras_L);
	// U
	frame_channel_u.FilteringOnNewImage(ras_U);
	// V
	frame_channel_v.FilteringOnNewImage(ras_V);
	// L+U
	frame_channel_lu.FilteringOnNewImage(ras_LU);
	// L+V
	frame_channel_lv.FilteringOnNewImage(ras_LV);
	// U+V
	frame_channel_uv.FilteringOnNewImage(ras_UV);
}

void MPLE_Luv::Synthesiz(Raster_Color<uchar> &image_output,const int kcluster)
{
	double				l,u,v;
	double				pLUV[3],prob;
	int					i,j,color_index;
	RF					rf_color_index;
	double				energy,energy_min;
	Raster_Color<uchar> image_color;
	RF3D				rf_iid_init;

	if (kcluster >= 0)
	{
		rf_iid_init.Set(frame_channel_l.GetLow(),frame_channel_l.GetHigh(),frame_channel_l.GetDelta(),
						frame_channel_u.GetLow(),frame_channel_u.GetHigh(),frame_channel_u.GetDelta(),
						frame_channel_v.GetLow(),frame_channel_v.GetHigh(),frame_channel_v.GetDelta());
		presp->ComputeHistogramAsIid(rf_iid_init,*p_mx_cluster_map,kcluster);
		rf_iid_init.SetCdf();
	}

	// set the present image to be a little bit bigger so that
	// boundary conditions can be avoided
	image_color.SetDimension(image_output.rows()+20,image_output.cols()+20);
	image_color.RGB.setTo(cv::Scalar::all(255));
	image_color.CreateLUVspace();

	// initilization
	rf_color_index.Set(0,num_total_color-1,1.0);
	for (i=0; i<num_total_color; i++)
		rf_color_index.SetOneBin(i,1.0);
	rf_color_index.SetCdf();

	for (j=0; j<image_color.rows(); j++)
		for (i=0; i<image_color.cols(); i++)
		{
			if (kcluster>=0)
			{
				rf_iid_init.Sampling(pLUV[0],pLUV[1],pLUV[2]);
				image_color.L(j,i) = (int)pLUV[0];
				image_color.U(j,i) = (int)pLUV[1];
				image_color.V(j,i) = (int)pLUV[2];
			}
			else
			{
				color_index = (int)rf_color_index.Sampling();
				image_color.L(j,i) = (int)p_important_color[0][color_index];
				image_color.U(j,i) = (int)p_important_color[1][color_index];
				image_color.V(j,i) = (int)p_important_color[2][color_index];
			}
		}

	int time=30,mm;
	FilteringOnNewImage(image_color);

	/* Jiadi
	if (pCtrProgressActive != NULL)
	{
		pCtrProgressActive->SetRange(0,time);
		pCtrProgressActive->SetPos(0);
		pCtrProgressActive->ShowWindow(SW_SHOW);
	}
	*/

	for (int t=0; t<time; t++)
	{
		/* Jiadi
		if (pCtrProgressActive != NULL)
		{
			pCtrProgressActive->SetPos(t);
			pCtrProgressActive->ShowWindow(SW_SHOW);
		}
		*/


		for (j=10; j<image_color.rows()-10; j++)
			for (i=10; i<image_color.cols()-10; i++)
			{
				// L
				frame_channel_l.ComputeEnergyAtEachValue(j,i);

				// U
				frame_channel_u.ComputeEnergyAtEachValue(j,i);
				// V
				frame_channel_v.ComputeEnergyAtEachValue(j,i);
				// L+U
				frame_channel_lu.ComputeEnergyAtEachValue(j,i);
				// L+V
				frame_channel_lv.ComputeEnergyAtEachValue(j,i);
				// U+V
				frame_channel_uv.ComputeEnergyAtEachValue(j,i);

				energy_min = 1e10;
				for (mm=0; mm<num_total_color; mm++)
				{
					l = p_important_color[0][mm];
					u = p_important_color[1][mm];
					v = p_important_color[2][mm];

					energy = Energy(l,u,v);
					prob = ProbFromEnergy(energy);
					if (energy<energy_min)
						energy_min = energy;

					rf_color_index.SetOneBin(mm,prob);
				}
				if (fabs(energy_min)>600)
				{
					for (mm=0; mm<num_total_color; mm++)
					{
						l = p_important_color[0][mm];
						u = p_important_color[1][mm];
						v = p_important_color[2][mm];

						energy = Energy(l,u,v)-energy_min;
						prob = ProbFromEnergy(energy);

						rf_color_index.SetOneBin(mm,prob);
					}
				}
					// sampling
				color_index = (int)rf_color_index.SamplingNonSetCdf();
				pLUV[0] = p_important_color[0][color_index];
				pLUV[1] = p_important_color[1][color_index];
				pLUV[2] = p_important_color[2][color_index];

				UpdateFilterRespAtOneSite(j,i,pLUV[0],pLUV[1],pLUV[2]);
			}

		for (j=0; j<image_color.rows(); j++)
			for (i=0; i<image_color.cols(); i++)
			{
				pLUV[0] = presp->filter_bank_L.ras_image(j,i);
				pLUV[1] = presp->filter_bank_U.ras_image(j,i);
				pLUV[2] = presp->filter_bank_V.ras_image(j,i);

				image_color.SetPixelInLUV(j,i,(int)pLUV[0],(int)pLUV[1],(int)pLUV[2]);
			}
		WriteToPPM(image_color,"aa.ppm");
	}

	for (j=0; j<image_output.rows(); j++)
		for (i=0; i<image_output.cols(); i++)
		{
			image_output.R(j,i) = image_color.R(j+10,i+10);
			image_output.G(j,i) = image_color.G(j+10,i+10);
			image_output.B(j,i) = image_color.B(j+10,i+10);
		}

	// hide progress window
	/* Jiadi
	if (pCtrProgressActive != NULL)
		pCtrProgressActive->ShowWindow(SW_HIDE);
	*/
}

// compute the probability map
void MPLE_Luv::ComputeEnergyMap(mcmcMatrix<double> &mx_prob_map)
{
	register int	i,j,m;
	double			l,u,v;
	double			energy,energy_min,Z;
	int				row_expansion=10,col_expansion=10;
	double			*penergys;

	// initilization
	mx_prob_map.SetDimension(presp->rows(),presp->cols());

	mx_prob_map.InitValue(0.0);
	presp->ExpandImage(row_expansion,col_expansion);
	presp->DoFiltering();

	// give more space than needed
	penergys = new double[num_total_color*2];

	/* Jiadi
	if (pCtrProgressActive != NULL)
	{
		pCtrProgressActive->SetRange(row_expansion,mx_prob_map.rows()+row_expansion);
		pCtrProgressActive->SetPos(0);
		pCtrProgressActive->ShowWindow(SW_SHOW);
	}
	*/

	// analytically compute the probobility
	for (j=row_expansion; j<mx_prob_map.rows()+row_expansion; j++)
		for (i=col_expansion; i<mx_prob_map.cols()+col_expansion; i++)
		{
			/* Jiadi
			if (pCtrProgressActive != NULL)
			{
				pCtrProgressActive->SetPos(j);
				pCtrProgressActive->ShowWindow(SW_SHOW);
			}
			*/

			// L
			frame_channel_l.ComputeEnergyAtEachValue(j,i);

			// U
			frame_channel_u.ComputeEnergyAtEachValue(j,i);
			// V
			frame_channel_v.ComputeEnergyAtEachValue(j,i);
			// L+U
			frame_channel_lu.ComputeEnergyAtEachValue(j,i);
			// L+V
			frame_channel_lv.ComputeEnergyAtEachValue(j,i);
			// U+V
			frame_channel_uv.ComputeEnergyAtEachValue(j,i);
				// compute the partition value
			energy_min = 1e10;
			for (m=0; m<num_total_color; m++)
			{
				l = p_important_color[0][m];
				u = p_important_color[1][m];
				v = p_important_color[2][m];

				penergys[m] = Energy(l,u,v);
				if (penergys[m]<energy_min)
					energy_min = penergys[m];
			}
			Z = 0.0;
			for (m=0; m<num_total_color; m++)
				Z += ProbFromEnergy(penergys[m]-energy_min);
			l = frame_channel_l.GetValueIndex(frame_channel_l.GetIndex(presp->ras_L(j-row_expansion,i-col_expansion)));
			u = frame_channel_u.GetValueIndex(frame_channel_u.GetIndex(presp->ras_U(j-row_expansion,i-col_expansion)));
			v = frame_channel_v.GetValueIndex(frame_channel_v.GetIndex(presp->ras_V(j-row_expansion,i-col_expansion)));
			energy = Energy(l,u,v);
			energy = energy-energy_min;

			mx_prob_map(j-row_expansion,i-col_expansion) = log(Z)+energy;
		}

	presp->ShrinkImageBack(row_expansion,col_expansion);
	presp->DoFiltering();

	if (penergys != NULL)	delete []penergys;
	/* Jiadi
	if (pCtrProgressActive != NULL)
		pCtrProgressActive->ShowWindow(SW_HIDE);
	*/
}


/////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////
MPLE_It::MPLE_It()
{
	bprecompute_resp = true;
	p_important_it_index = NULL;
	p_important_it = NULL;
}

MPLE_It::~MPLE_It()
{
	if (p_important_it_index != NULL)
		delete [](p_important_it_index);
	if (p_important_it != NULL)
		delete [](p_important_it);
}

bool MPLE_It::ValidIt(const double it)
{
	int k_it;

	k_it  = frame_channel_it.GetIndex(it);

	if (frame_channel_it.ValidIndex(k_it))
		return true;
	else
		return false;
}

void MPLE_It::Initialization(void)
{
	frame_channel_it.Set(presp->filter_bank_intensity,bprecompute_resp,biggest_rowsize,biggest_colsize);

	int			j,m;
	double		it;
	Histogram	*phist_it;

	phist_it = &frame_channel_it.histogram_ref_input;
	num_total_it = phist_it->count;

	m = 0;
	if (p_important_it_index != NULL)
		delete []p_important_it_index;
	p_important_it_index = new int[num_total_it];
	if (p_important_it != NULL)
		delete []p_important_it;
	p_important_it = new double[num_total_it];

	for (j=0; j<phist_it->count; j++)
	{
		it = phist_it->GetValueIndex(j);

		p_important_it[m] = it;
		p_important_it_index[m] = frame_channel_it.GetIndex(it);

		m++;
	}
}

bool MPLE_It::UpdateBeta(const double d)
{
	// L
	return frame_channel_it.UpdateBeta(d);
}

void MPLE_It::ComputeSiteHistogramInAdvance(const int kcluster)
{
	if (bprecompute_resp)
	{
		frame_channel_it.ComputeResponse(presp->cut_off,presp->rows()-(presp->cut_off+1),
										presp->cut_off,presp->cols()-(presp->cut_off+1));

//		SaveRespHistToFile(0);
	}
}

void MPLE_It::SaveRespHistToFile(const int m)
{
	char		str[200];
	ofstream	file;

	sprintf(str,"d:\\temp\\%d.rsp", m);
	file.open(str, ios::out| ios::binary);
	if (file.fail())
		return;

	frame_channel_it.SaveRespHistToFile(file);

	file.close();
}

void MPLE_It::ReadRespHistFromFile(const int m)
{
	if (current_patch != m)
	{
		char		str[200];
		ifstream	file;

		sprintf(str,"d:\\temp\\%d.rsp", m);

		file.open(str, ios::in| ios::binary|ios::_Nocreate);
		if (file.fail())
			return;

		frame_channel_it.ReadRespHistFromFile(file);

		file.close();

		current_patch = m;
	}
}

double	MPLE_It::Energy(const double it)
{
	int			k;
	double		energy;

	energy = 0.0;
	// L
	k = frame_channel_it.GetIndex(it);
	if (frame_channel_it.vc_energy.Valid(k))
		energy += frame_channel_it.vc_energy(k);
	else
		return 9000.0;

	return energy;
}

double MPLE_It::Energy(const int k_it)
{
	double		energy;

	energy = frame_channel_it.vc_energy(k_it);

	return energy;
}


void MPLE_It::UpdateHistogramSyn(void)
{
	frame_channel_it.UpdateHistogramSyn();
}

// from precomputed histogram
void MPLE_It::UpdateHistogramSynTemp(const double prob, const double it, const int row, const int col)
{
	// L
	frame_channel_it.UpdateHistogramSynTemp(prob,it,row,col);
}

void MPLE_It::ClearHistogramSyn(void)
{
	// L
	frame_channel_it.ClearHistogramSyn();
}

void MPLE_It::ClearHistogramSynTemp(void)
{
	// L
	frame_channel_it.ClearHistogramSynTemp();
}

void MPLE_It::UpdateFilterRespAtOneSite(const int j, const int i, const double it)
{
	// L
	frame_channel_it.UpdateFilterRespAtOneSite(j,i,it);
}


void MPLE_It::Learn(const int time, const bool bprecompute, const int kcluster)
{
	register int	i,j,t;
	register double	it;
	double			prob,energy,energy_min;
	mcmcVector<double>	vc_energy;
	int				row_size,col_size;

	// precompute all the responses regarding to the changes of pixel values
	bprecompute_resp = bprecompute;

	Initialization();
	// temporary energy
	vc_energy.SetRowVector(num_total_it);

	current_patch = -1;

	char btafile_name[200];
	sprintf(btafile_name,"a_%d.bta",kcluster);

	if (time > 0)
	{
		ComputeSiteHistogramInAdvance(kcluster);

		/* Jiadi
		if (pCtrProgressActive != NULL)
		{
			pCtrProgressActive->SetRange(0,time);
			pCtrProgressActive->SetPos(0);
			pCtrProgressActive->ShowWindow(SW_SHOW);
		}
		*/
	}

	row_size = presp->rows()-2*presp->cut_off;
	col_size = presp->cols()-2*presp->cut_off;

	register int mm;
	bool	bConverge=false;

	t = 0;
	while (t<time && !bConverge)
	{
		ClearHistogramSyn();

		for (j=0; j<row_size; j++)
		{
			for (i=0; i<col_size; i++)
			{
				if (p_mx_cluster_map->Data(j+presp->cut_off,i+presp->cut_off)==kcluster)
				{
					ClearHistogramSynTemp();
					energy_min = 1e10;

					// L
					frame_channel_it.SaveTempEnergy(j,i);	// from the pre-computed histog

					// get the minimal energy value
					for (mm=0; mm<num_total_it; mm++)
					{
						energy = Energy(p_important_it_index[mm]);
						vc_energy(mm) = energy;
							if (energy<energy_min)	energy_min = energy;
					}
					for (mm=0; mm<num_total_it; mm++)
					{
						it = p_important_it[mm];
							energy = vc_energy(mm)-energy_min;
						prob = ProbFromEnergy(energy);
						UpdateHistogramSynTemp(prob,it,j,i);
					}
					UpdateHistogramSyn();
				}
			}
		}
		// update betas
		bConverge = UpdateBeta(GetD(t));
		t++;
		if (t>=time)
			break;
		/* Jiadi
		if (pCtrProgressActive != NULL)
		{
			// show progress
			pCtrProgressActive->SetPos(t);
			pCtrProgressActive->ShowWindow(SW_SHOW);
		}
		*/
		if (bConverge)
			break;
		SaveBetaToFile(btafile_name);
		if (bConverge)
			break;
	}
	/* Jiadi
	if (pCtrProgressActive != NULL)
		pCtrProgressActive->ShowWindow(SW_HIDE);
	*/

	if (time != 0)
		SaveBetaToFile(btafile_name);
	ReadBetaFromFile(btafile_name);
}

void MPLE_It::SaveBetaToFile(char *file_name)
{
	FILE		*pfile;
	Histogram	*phist;
	int			i;

	if (file_name == NULL)
		return;

	pfile = fopen(file_name,"w+");
	if(pfile != NULL)
	{
		for (i=0; i<frame_channel_it.vc_B.size(); i++)
		{
			phist = (Histogram *)&frame_channel_it.vc_B(i);
			phist->WriteToTxtFile(*pfile);
		}
		fprintf(pfile,"-----------------------\n");
		fclose(pfile);
	}
}

void MPLE_It::ReadBetaFromFile(char *file_name)
{
	FILE		*pfile;
	Histogram	*phist;
	int			i;
	char		str[100];
	if (file_name == NULL)
		return;

	pfile = fopen(file_name,"r");
	if (pfile != NULL)
	{
		for (i=0; i<frame_channel_it.vc_B.size(); i++)
		{
			phist = (Histogram *)&frame_channel_it.vc_B(i);
			phist->ReadFromTxtFile(*pfile);
		}
		fscanf(pfile,"%s",&str);
		fclose(pfile);
	}
}

void MPLE_It::FilteringOnNewImage(Raster<uchar> &image_intensity)
{
	Raster<double> ras_intensity;

	ras_intensity.SetDimension(image_intensity.rows(),image_intensity.cols());

	for (int j=0; j<ras_intensity.rows(); j++)
		for (int i=0; i<ras_intensity.cols(); i++)
			ras_intensity(j,i) = (double)image_intensity(j,i);

	// L
	frame_channel_it.FilteringOnNewImage(ras_intensity);
}

void MPLE_It::Synthesiz(Raster<uchar> &image_output,const int kcluster)
{
	double				it;
	double				prob;
	register int		i,j,it_index;
	RF					rf,rf_it_index;
	double				energy,energy_min;
	Raster<uchar>		image;
	RF					rf_iid_init;

	if (kcluster >= 0)
	{
		rf_iid_init.Set(frame_channel_it.GetLow(),frame_channel_it.GetHigh(),frame_channel_it.GetDelta());
		presp->ComputeHistogramAsIid(rf_iid_init,*p_mx_cluster_map,kcluster);
		rf_iid_init.SetCdf();
	}

	// set the present image to be a little bit bigger so that
	// boundary conditions can be avoided
	image.SetDimension(image_output.rows()+20,image_output.cols()+20);
	// initilization
	rf_it_index.Set(0,num_total_it-1,1.0);
	for (i=0; i<num_total_it; i++)
		rf_it_index.SetOneBin(i,1.0);
	rf_it_index.SetCdf();

	rf.Set(frame_channel_it.GetLow(),frame_channel_it.GetHigh(),frame_channel_it.GetDelta());

	for (j=0; j<image.rows(); j++)
		for (i=0; i<image.cols(); i++)
		{
			if (kcluster>=0)
				image(j,i) = (uchar)rf_iid_init.Sampling();
			else
			{
				it_index = (int)rf_it_index.Sampling();
				image(j,i) = (uchar)p_important_it[it_index];
			}
		}

	int time=30,mm;
	FilteringOnNewImage(image);

	/* Jiadi
	if (pCtrProgressActive != NULL)
	{
		pCtrProgressActive->SetRange(0,time);
		pCtrProgressActive->SetPos(0);
		pCtrProgressActive->ShowWindow(SW_SHOW);
	}
	*/

	for (int t=0; t<time; t++)
	{
		/* Jiadi
		if (pCtrProgressActive != NULL)
		{
			pCtrProgressActive->SetPos(t);
			pCtrProgressActive->ShowWindow(SW_SHOW);
		}
		*/

		for (j=10; j<image.rows()-10; j++)
			for (i=10; i<image.cols()-10; i++)
			{
				// L
				frame_channel_it.ComputeEnergyAtEachValue(j,i);

				energy_min = 1e10;
				for (mm=0; mm<num_total_it; mm++)
				{
					energy = Energy(p_important_it_index[mm]);
					if (energy<energy_min)
						energy_min = energy;
				}
				for (mm=0; mm<num_total_it; mm++)
				{
					energy = Energy(p_important_it_index[mm])-energy_min;
					prob = ProbFromEnergy(energy);

					rf_it_index.SetOneBin(mm,prob);
				}
				// sampling
				it_index = (int)rf_it_index.SamplingNonSetCdf();
				if (it_index>0 && it_index<num_total_it)
				{
					it = p_important_it[it_index];
					UpdateFilterRespAtOneSite(j,i,it);
				}
			}
		for (j=0; j<image.rows(); j++)
			for (i=0; i<image.cols(); i++)
				image(j,i) = (uchar)presp->filter_bank_intensity.ras_image(j,i);

		WriteToPPM(image,"aa.ppm");
	}

	for (j=0; j<image_output.rows(); j++)
		for (i=0; i<image_output.cols(); i++)
			image_output(j,i) = image(j+10,i+10);

	// hide progress window
	/* Jiadi
	if (pCtrProgressActive != NULL)
		pCtrProgressActive->ShowWindow(SW_HIDE);
	*/
}


void MPLE_It::ComputeEnergyMap(mcmcMatrix<double> &mx_prob_map)
{
	register int	i,j;
	double			energy,energy_min,Z;
	int				row_expansion=10,col_expansion=10;
	double			*penergys;

	double			it;
	int				mm;
	Vector_Double	vc_resp;

	// initilization
	mx_prob_map.SetDimension(presp->rows(),presp->cols());
	mx_prob_map.InitValue(0.0);
	presp->ExpandImage(row_expansion,col_expansion);
	presp->DoFiltering();

	// give more space than needed
	penergys = new double[num_total_it];

	/* Jiadi
	if (pCtrProgressActive != NULL)
	{
		pCtrProgressActive->SetRange(row_expansion,mx_prob_map.rows()+row_expansion);
		pCtrProgressActive->SetPos(0);
		pCtrProgressActive->ShowWindow(SW_SHOW);
	}
	*/

	// analytically compute the probobility
	for (j=row_expansion; j<mx_prob_map.rows()+row_expansion; j++)
		for (i=col_expansion; i<mx_prob_map.cols()+col_expansion; i++)
		{
			/* Jiadi
			if (pCtrProgressActive != NULL)
			{
				pCtrProgressActive->SetPos(j);
				pCtrProgressActive->ShowWindow(SW_SHOW);
			}
			*/

			// compute the partition value
			energy_min = 1e10;

			frame_channel_it.ComputeEnergyAtEachValue(j,i);

			for (mm=0; mm<num_total_it; mm++)
			{
				it = p_important_it[mm];
				penergys[mm] = Energy(it);

				if (penergys[mm]<energy_min)
					energy_min = penergys[mm];
			}

			Z = 0.0;
			for (mm=0; mm<num_total_it; mm++)
				Z += ProbFromEnergy(penergys[mm]-energy_min);
			energy = Energy(presp->ras_intensity(j-row_expansion,i-col_expansion))-energy_min;

			mx_prob_map(j-row_expansion,i-col_expansion) = log(Z)+energy;
		}

	presp->ShrinkImageBack(row_expansion,col_expansion);
	presp->DoFiltering();

	if (penergys != NULL)	delete []penergys;
	/* Jiadi
	if (pCtrProgressActive != NULL)
		pCtrProgressActive->ShowWindow(SW_HIDE);
	*/
}

/////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////
// learn all the pseudo likelihood parameters for a given clustered image

// construction & destruction
MCMC_Learn_PLE_Clusters::MCMC_Learn_PLE_Clusters()
{
	num_clusters		= 0;
	p_cluster_center	= NULL;
	pbound				= NULL;
}

MCMC_Learn_PLE_Clusters::~MCMC_Learn_PLE_Clusters()
{
	if (p_cluster_center != NULL)
		delete []p_cluster_center;
	p_cluster_center = NULL;

	if (pbound != NULL)
		delete []pbound;
	pbound = NULL;
}

// operations
mcmcMatrix<int>	MCMC_Learn_PLE_Clusters::GetOptimalDivision(const int size_whole, const int size_patch,
														const int size_bound_cut, const int size_overlap,
														int *p_size_biggest)
{
	int			num,size_left,current_pos;
	mcmcMatrix<int>	mx_par_temp,mx_partition;

	num = (int)ceil(((float)size_whole)/size_patch);
	mx_par_temp.SetDimension(num,2);

	num = 0;
	size_left = size_whole-2*size_bound_cut;
	current_pos = size_bound_cut;
	while (size_left>0)
	{
		if (size_left>size_patch*3/2)
		{
			// starting position of this patch
			mx_par_temp(num,0) = current_pos;
			// ending position of this patch
			current_pos += size_patch;
			mx_par_temp(num,1) = size_patch+size_overlap;
			size_left -= size_patch;
		}
		else
		{
			// starting position of this patch
			mx_par_temp(num,0) = current_pos;
			mx_par_temp(num,1) = size_left;
			size_left = 0;
		}
		num++;
	}

	int size_biggest=0;
	mx_partition.SetDimension(num,2);
	for (int i=0; i<num; i++)
	{
		mx_partition(i,0) = mx_par_temp(i,0);
		mx_partition(i,1) = mx_par_temp(i,1);

		if (size_biggest<mx_partition(i,1))
			size_biggest = mx_partition(i,1);
	}
	if (p_size_biggest != NULL)
		*p_size_biggest = size_biggest;

	return mx_partition;
}


void MCMC_Learn_PLE_Clusters::Set(mcmcMatrix<int> &mx_map, MCMC_FilterResp &resp1)
{
	int i,j;

	presp = &resp1;
	mx_cluster_map = mx_map;
	// how many clusters are there
	num_clusters = mx_cluster_map.GetMax()+1;
	// initialization
	vc_valid_clusters.SetDimension(num_clusters,1);
	vc_size_clusters.SetRowVector(num_clusters);
	vc_size_clusters.InitValue(0);

	// number of pixles for each cluster
	for (j=0; j<mx_cluster_map.rows(); j++)
		for (i=0; i<mx_cluster_map.cols(); i++)
			vc_size_clusters(mx_cluster_map(j,i))++;
	for (i=0; i<vc_size_clusters.size(); i++)
	{
		if (vc_size_clusters(i)<60*60)
			vc_valid_clusters(i) = false;
		else
			vc_valid_clusters(i) = true;
	}
}


void MCMC_Learn_PLE_Clusters::LearnPLE4Clusters(Gibbs_B *pref_l, Gibbs_B *pref_u, Gibbs_B *pref_v,
												mcmcMatrix<double> *pmx_energy_map, Raster_Color<uchar> *pras_syn)
{
	Raster_Color<uchar> image_syn;
	int i,j,m,n,k;

	biggest_rowsize = 100;
	biggest_colsize = 100;

	// set pseudo
	mple_luv.Set(*presp,mx_cluster_map,mx_partition,biggest_rowsize,biggest_colsize);

	PrepareForPasting();

	for (k=0; k<num_clusters; k++)
		if (vc_valid_clusters(k))
		{

			PasteAnImageToLearn(k);
			presp->ComputeHistogram4Mask(mx_cluster_map, k);
			mple_luv.Learn(500,true,k);

			pref_l[k] = mple_luv.frame_channel_l.vc_B(0);
			pref_u[k] = mple_luv.frame_channel_u.vc_B(0);
			pref_v[k] = mple_luv.frame_channel_v.vc_B(0);

			// compute energies for each pixel
			if (pmx_energy_map != NULL)
				mple_luv.ComputeEnergyMap(pmx_energy_map[k]);

			char name[200];
			sprintf(name,"aba_%d.eng",k);
//			if (pmx_energy_map != NULL)
//				WriteBigMatrix(pmx_energy_map[k],name);

			// synthesize a color image
			if (pras_syn != NULL)
			{
				image_syn.SetDimension(50,50);
				mple_luv.Synthesiz(image_syn,k);
				// paste synthesized image
				pras_syn[k].SetDimension(presp->rows(),presp->cols());
				for (j=0; j<pras_syn[k].rows(); j++)
					for (i=0; i<pras_syn[k].cols(); i++)
					{
						m = 1+j%48;
						n = 1+i%48;
						pras_syn[k].R(j,i) = image_syn.R(m,n);
						pras_syn[k].G(j,i) = image_syn.G(m,n);
						pras_syn[k].B(j,i) = image_syn.B(m,n);
					}
			}

			sprintf(name,"aba_syn_%d.ppm",k);
//			WriteToPPM(pras_syn[k],name);

			// the original image has been destroyed by the synthesis process
			// copy it back
			mple_luv.FilteringOnNewImage(presp->image_original);
		}
}

void MCMC_Learn_PLE_Clusters::PrepareForPasting(void)
{
	int i,j,n;

	// reset labels of the clustering map to have connected clusters
	ras_connected_cluster = mx_cluster_map;
	num_connected_cluster = ras_connected_cluster.Relabel(false)+1;
	if (p_cluster_center != NULL)
		delete []p_cluster_center;
	p_cluster_center = new McPoint2D[num_connected_cluster];
	vc_size_connected.SetRowVector(num_connected_cluster);
	vc_cluster_index.SetRowVector(num_connected_cluster);
	if (pbound != NULL)
		delete []pbound;
	pbound = new irmm_Rect[num_connected_cluster];

	// initialization
	for (i=0; i<num_connected_cluster; i++)
	{
		vc_size_connected(i) = 0;
		p_cluster_center[i].Set(0,0);
		pbound[i].Set(9999,9999,-1,-1);
	}

	for (j=0; j<mx_cluster_map.rows(); j++)
		for (i=0; i<mx_cluster_map.cols(); i++)
		{
			// which connected part
			n = ras_connected_cluster(j,i);

			vc_size_connected(n)++;
			p_cluster_center[n].x += j;
			p_cluster_center[n].y += i;
			vc_cluster_index(n) = mx_cluster_map(j,i);
			pbound[n].SetBound(i,j);

		}
	for (n=0; n<num_connected_cluster; n++)
	{
		p_cluster_center[n].x = (int)(p_cluster_center[n].x/vc_size_connected(n));
		p_cluster_center[n].y = (int)(p_cluster_center[n].y/vc_size_connected(n));
	}

	rf_connected_cluster.Set(0,num_connected_cluster-1,1);
}

void MCMC_Learn_PLE_Clusters::PasteAnImageToLearn(const int kcluster)
{
	int				i,j,n;
	RF3D			rf_pixel;
	mcmcMatrix<char>	mask_pixel;

	for (n=0; n<num_connected_cluster; n++)
	{
		if (vc_cluster_index(n)==kcluster)
			rf_connected_cluster.SetOneBin(n,vc_size_connected(n));
		else
			rf_connected_cluster.SetOneBin(n,0);
	}

	// pixels that are unfilled
	mask_pixel.SetDimension(biggest_rowsize,biggest_colsize);
	mask_pixel.InitValue(0);
	image_paste.SetDimension(mask_pixel.rows(),mask_pixel.cols());
	rf_pixel.Set(0,mask_pixel.rows()-1,1,
		         0,mask_pixel.cols()-1,1,
				 0,1,1);
	for (j=0; j<mask_pixel.rows(); j++)
		for (i=0; i<mask_pixel.cols(); i++)
		{
			rf_pixel.SetOneBin(j,i,0,1.0);
			rf_pixel.SetOneBin(j,i,1,1.0);
		}

	double	drow,dcol,dump;
	int		rowc,colc,row,col;
	int	filled=0;

	while (filled < mask_pixel.size())
	{
		// pick up one patch
		int patch=(int)rf_connected_cluster.SamplingNonSetCdf();
		// pick up one pixel
		rf_pixel.SamplingNonSetCdf(drow,dcol,dump);
		rowc = (int)drow; colc = (int)dcol;

		for (j=pbound[patch].top; j<=pbound[patch].bottom; j++)
			for (i=pbound[patch].left; i<=pbound[patch].right; i++)
			{
				row = j-(int)(p_cluster_center[patch].x)+rowc;
				col = i-(int)(p_cluster_center[patch].y)+colc;

				if (mask_pixel.Valid(row,col) && mask_pixel(row,col)==0 && mx_cluster_map(j,i)==kcluster)
				{
					mask_pixel(row,col) = 1;
					image_paste.R(row,col) = presp->image_original.R(j,i);
					image_paste.G(row,col) = presp->image_original.G(j,i);
					image_paste.B(row,col) = presp->image_original.B(j,i);

					rf_pixel.SetOneBin(row,col,0,0);
					rf_pixel.SetOneBin(row,col,1,0);

					filled++;
				}
			}
	}

	char name[200];

	sprintf(name,"c_%d.ppm",kcluster);

	WriteToPPM(image_paste,name);
}


void MCMC_Learn_PLE_Clusters::LearnPLE4Clusters(Gibbs_B *pref_it,mcmcMatrix<double> *pmx_energy_map, Raster<uchar> *pras_syn)
{
	Raster<uchar>	image_syn;
	int				i,j,m,n,k;
	bool			bComputeSiteResp=true;

	biggest_rowsize = presp->rows();
	biggest_colsize = presp->cols();
	// set pseudo
	mple_it.Set(*presp,mx_cluster_map,mx_partition,biggest_rowsize,biggest_colsize);

	for (k=0; k<num_clusters; k++)
		if (vc_valid_clusters(k))
		{
			// Jiadi Yang -----------------------------
			Raster<uchar> rRGB;
			cv::Mat mRGB[3];
			cv::split(presp->image_original.RGB, mRGB);
			rRGB.SetDimension(presp->image_original.rows(), presp->image_original.cols());
			rRGB.data = mRGB[0];
			
			presp->FilteringOnNewImage(rRGB);
			presp->ComputeHistogram4Mask(mx_cluster_map,k);

			// ----------------------------------------
			mple_it.Learn(100,bComputeSiteResp,k);

			// share the same response values
			bComputeSiteResp = false;

			pref_it[k] = mple_it.frame_channel_it.vc_B(0);

			// compute energies for each pixel
			if (pmx_energy_map != NULL)
				mple_it.ComputeEnergyMap(pmx_energy_map[k]);

			// synthesize a color image
			if (pras_syn != NULL)
			{
				image_syn.SetDimension(50,50);
				mple_it.Synthesiz(image_syn,k);
				// paste synthesized image
				pras_syn[k].SetDimension(presp->image_original.rows(),presp->image_original.cols());
				for (j=0; j<pras_syn[k].rows(); j++)
					for (i=0; i<pras_syn[k].cols(); i++)
					{
						m = 1+j%48;
						n = 1+i%48;
						pras_syn[k](j,i) = image_syn(m,n);
					}
			}
			// the original image has been destroyed by the synthesis process
			// copy it back
			mple_it.FilteringOnNewImage(rRGB);
		}
}

bool WriteToPPM(Raster_Color<uchar> &image_color, char *filename)
{
	FILE *file;
	file = fopen(filename,"w+");
	if (file != NULL)
	{
		fprintf(file,"P3 %d %d 255\n", image_color.cols(),image_color.rows());
		for (int j=0; j<image_color.rows(); j++)
		{
			for (int i=0; i<image_color.cols(); i++)
				fprintf(file, "%d %d %d ", image_color.R(j,i),image_color.G(j,i),image_color.B(j,i));
			fprintf(file,"\n");
		}
		fclose(file);
		return true;
	}
	else
		return false;
}

bool WriteToPPM(Raster<uchar> &image, char *filename)
{
	FILE *file;
	file = fopen(filename,"w+");
	if (file != NULL)
	{
		fprintf(file,"P3 %d %d 255\n", image.cols(),image.rows());
		for (int j=0; j<image.rows(); j++)
		{
			for (int i=0; i<image.cols(); i++)
				fprintf(file, "%d %d %d ", image(j,i),image(j,i),image(j,i));
			fprintf(file,"\n");
		}
		fclose(file);
		return true;
	}
	else
		return false;
}

bool ReadFromPPM(Raster_Color<uchar> &image_color, char *file_name)
{
	FILE *file;
	int nrow,ncol;
	int r,g,b;

	file = fopen(file_name,"r");
	if (file != NULL)
	{
		fscanf(file,"P3 %d %d 255", &ncol,&nrow);
		image_color.SetDimension(nrow,ncol);
		for (int j=0; j<image_color.rows(); j++)
		{
			for (int i=0; i<image_color.cols(); i++)
			{
				fscanf(file, "%d %d %d", &r,&g, &b);
				image_color.R(j,i) = r;
				image_color.G(j,i) = g;
				image_color.B(j,i) = b;
			}
		}
		fclose(file);
		return true;
	}
	else
		return false;
}

bool ReadFromPPM(Raster<uchar> &image, char *file_name)
{
	FILE	*file;
	int		nrow,ncol,dummy1,dummy2,value;

	file = fopen(file_name,"r");
	if (file != NULL)
	{
		fscanf(file,"P3 %d %d 255", &ncol,&nrow);
		image.SetDimension(nrow,ncol);
		for (int j=0; j<image.rows(); j++)
		{
			for (int i=0; i<image.cols(); i++)
			{
				fscanf(file, "%d %d %d ", &value,&dummy1,&dummy2);
				image(j,i) = (uchar)value;
			}
		}
		fclose(file);
		return true;
	}
	else
		return false;
}
